#include "SsimCompare.h"

cv::Scalar SsimCompare::matCompare(std::string matSrcPath, std::string matDstPath)
{
    cv::Mat matSrc = cv::imread(matSrcPath);
    cv::Mat matDst = cv::imread(matDstPath);

    if (matSrc.empty() || matDst.empty())
    {
        return 0.0;
    }

    cv::resize(matDst, matDst, cv::Size(matSrc.cols, matSrc.rows));

    auto start = std::chrono::steady_clock::now();

	cv::Scalar result = this->getMSSIM(matSrc, matDst);

    auto end = std::chrono::steady_clock::now();

    std::cout << "Time of PSNR CPU " << std::chrono::duration<double, std::micro>(end - start).count() << " milliseconds."
        << " With result of B" << result.val[0] << " G" << result.val[1] << " R" << result.val[2] << std::endl; ;

    return result;
}

cv::Scalar SsimCompare::getMSSIM(const cv::Mat& matSrc, const cv::Mat& matDst)
{
	const double C1 = 6.5025, C2 = 58.5225;
	/***************************** INITS **********************************/
	int d = CV_32F;

	cv::Mat I1, I2;
	matSrc.convertTo(I1, d);           // cannot calculate on one byte large values
	matDst.convertTo(I2, d);

	cv::Mat I2_2 = I2.mul(I2);        // I2^2
	cv::Mat I1_2 = I1.mul(I1);        // I1^2
	cv::Mat I1_I2 = I1.mul(I2);        // I1 * I2

	/*************************** END INITS **********************************/

	cv::Mat mu1, mu2;   // PRELIMINARY COMPUTING
	cv::GaussianBlur(I1, mu1, cv::Size(11, 11), 1.5);
	cv::GaussianBlur(I2, mu2, cv::Size(11, 11), 1.5);

	cv::Mat mu1_2 = mu1.mul(mu1);
	cv::Mat mu2_2 = mu2.mul(mu2);
	cv::Mat mu1_mu2 = mu1.mul(mu2);

	cv::Mat sigma1_2, sigma2_2, sigma12;

	cv::GaussianBlur(I1_2, sigma1_2, cv::Size(11, 11), 1.5);
	sigma1_2 -= mu1_2;

	cv::GaussianBlur(I2_2, sigma2_2, cv::Size(11, 11), 1.5);
	sigma2_2 -= mu2_2;

	cv::GaussianBlur(I1_I2, sigma12, cv::Size(11, 11), 1.5);
	sigma12 -= mu1_mu2;

	// FORMULA
	cv::Mat t1, t2, t3;

	t1 = 2 * mu1_mu2 + C1;
	t2 = 2 * sigma12 + C2;
	t3 = t1.mul(t2);              // t3 = ((2*mu1_mu2 + C1).*(2*sigma12 + C2))

	t1 = mu1_2 + mu2_2 + C1;
	t2 = sigma1_2 + sigma2_2 + C2;
	t1 = t1.mul(t2);               // t1 =((mu1_2 + mu2_2 + C1).*(sigma1_2 + sigma2_2 + C2))

	cv::Mat ssim_map;
	cv::divide(t3, t1, ssim_map);      // ssim_map =  t3./t1;

	cv::Scalar mssim = cv::mean(ssim_map); // mssim = average of ssim map

	return mssim;
}
